Post-unblinding augmentation trial data

library(tidyverse)
library(magrittr)
library(gt)
library(patchwork)
library(SummarizedExperiment)
library(tidySummarizedExperiment)
library(MultiAssayExperiment)

tmp <- fs::dir_map("R scripts/", source)
tmp <- fs::dir_map("../VIBRANT-99-utils/R/", source)
rm(tmp)
theme_set(theme_light())

The purpose of this document is to compute additional “calculated variables” that require unblinded data, and to save the updated MultiAssayExperiment object for downstream analyses.

Loading the MultiAssayExperiment object

mae_files <- 
  fs::dir_ls(
    get_04_output_dir(),
    regexp = "mae_full_.*\\.rds$"
  ) |> 
  sort(decreasing = TRUE) |>
  magrittr::extract(1)

mae <- readRDS(mae_files)

rm(mae_files)

Calculated variables

Study product exposure related variables

Exposure by qPCR on dosing days

After Unblinding, qPCR assays were performed on daily self-collected swabs to inform effective exposure (see 30_qPCR_daily.qmd).

In this section, this daily qPCR data is analysed to describe patterns of LBP detection across the dosing period and subsequent follow-up.

df <- 
  mae[["qPCR_daily"]] |> 
  as_tibble() |> 
  select(uid, .feature, copies_per_swab_med, taxon_label, LBP, strain_origin) |> 
  left_join(
    mae@colData |> 
      as_tibble() |> 
      select(
        uid, 
        pid, 
        visit_code, PIPV, visit, study_day, 
        site, randomized, arm, ITT, mITT, PP
      ),
    by = join_by(uid)
  ) |> 
  arrange(LBP, strain_origin) |>
  mutate(.feature = .feature |> fct_inorder())

16S total copies per swabs

First, to assess general quality of collected swabs, we visualize the qPCR 16S results across visits and participants. For each participant, the median number of 16S copies per swab (log₁₀-transformed) is displayed as a heatmap, with participants on the y-axis and visit codes on the x-axis. The plot is stratified by study arm and site using nested facets.

df |> 
  filter(.feature == "16S", randomized) |> 
  ggplot() +
  aes(x = visit_code, y = pid, fill = copies_per_swab_med |> log10()) +
  ggh4x::facet_nested(arm + site ~ ., scales = "free", space = "free") +
  geom_tile() +
  scale_fill_viridis_c(option = "F", na.value = "blue") +
  theme(
    axis.text = element_text(size = 6),
    axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 6)
  )

LBP strains

We generate heatmaps of qPCR results for each study arm. For every participant and visit, the median number of copies per swab for each LBP strain is shown, with strain identity on the y-axis and visit code on the x-axis.

plots <- 
  map(
  .x = df$arm |> unique() |> sort(),
  ~ df |> 
    filter(.feature != "16S", arm == .x) |> 
    ggplot() +
    aes(x = visit_code, y = .feature, fill = LBP, alpha = copies_per_swab_med |> asinh()) +
    ggh4x::facet_nested(site + pid ~ ., scales = "free", space = "free") +
    geom_tile() +
    scale_alpha_continuous(range = c(0.1, 1), limits = c(-1/2, 20)) +
    ggtitle(.x) +
    theme(
      axis.text = element_text(size = 4),
      axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 6)
    )
)

names(plots) <- df$arm |> unique() |> sort()

Pl

LC106-7

LC106-3

LC106-o

LC115

We plot longitudinal qPCR measurements for each participant and strain, stratified by study arm, site, and participant. The x-axis shows study days and the y-axis the median copies per swab (log₁₀ scale). Each line represents one LBP strain, with color indicating strain identity. Shaded background areas indicate treatment periods: red for oral metronidazole and yellow for LBP dosing. A dashed horizontal line marks the 10⁵ copies/swab threshold.

plots <- 
  map(
  .x = df$arm |> unique() |> sort(),
  ~ df |> 
    filter(.feature != "16S", arm == .x, study_day >= 0) |> 
    ggplot() +
    aes(x = study_day, y = copies_per_swab_med, col = LBP ) + #  |> asinh()
    ggh4x::facet_nested(site + pid ~ .) + # scales = "free", space = "free"
    geom_hline(yintercept = 10^5, lty = 2, col = "gray40") +
    geom_rect(
      data = tibble(study_day = 0, copies_per_swab_med = NA), 
      xmin = 0.5, xmax = 7.5, ymin = -Inf, ymax = Inf, alpha = 0.1, fill = "red", col = NA
    ) +
    geom_rect(
      data = tibble(study_day = 0, copies_per_swab_med = NA), 
      xmin = 7.5, 
      xmax = case_when(
        .x == "LC106-3" ~ 10.5,
        .x == "LC106-o" ~ 10.5,
        TRUE ~ 14.5
      ), 
      ymin = -Inf, ymax = Inf, alpha = 0.1, fill = "gold", col = NA
    ) +
    # Rectangle jaune supplémentaire pour LC106-o uniquement
    {if(.x == "LC106-o") geom_rect(
      data = tibble(study_day = 0, copies_per_swab_med = NA), 
      xmin = 3.5, xmax = 7.5, ymin = -Inf, ymax = Inf, alpha = 0.1, fill = "gold", col = NA
    )} +
    geom_line(aes(group = .feature)) +
    geom_point() +
    scale_y_log10() +
    expand_limits(y = 10^10) +
    scale_x_continuous(breaks = seq(0, 35, by = 7)) +
    ggtitle(.x)
)

names(plots) <- df$arm |> unique() |> sort()

Pl

LC106-7

LC106-3

LC106-o

LC115

LBP exposure

We summarize the qPCR results to compute total LBP copies per swab and the number of detected LBP strains for each participant at each visit. We also calculate the proportion of detected LBP strains among those expected based on the study arm.

df_summ <- 
  df |> 
  # filter(.feature != "16S") |> 
  group_by(uid, pid, randomized, arm, ITT, mITT, PP, site, visit_code, PIPV, visit, study_day) |> 
  summarise(
    LC106_total = sum(copies_per_swab_med[str_detect(LBP, "LC106")], na.rm = TRUE),
    LC115_excl_total = sum(copies_per_swab_med[!str_detect(LBP, "LC106")], na.rm = TRUE),
    total_LBP = sum(copies_per_swab_med[str_detect(LBP, "LC")], na.rm = TRUE),
    total_16S = sum(copies_per_swab_med[.feature == "16S"], na.rm = TRUE),
    n_LC106 = sum(str_detect(LBP, "LC106") & (copies_per_swab_med != 0), na.rm = TRUE),
    n_LC115_excl = sum(!str_detect(LBP, "LC106") & (copies_per_swab_med != 0), na.rm = TRUE),
    n_LBP = sum(str_detect(LBP, "LC") & (copies_per_swab_med != 0), na.rm = TRUE),
    .groups = "drop"
  ) |> 
  mutate(
    n_expected = 
      case_when(
        arm == "Pl" ~ 0,
        str_detect(arm, "LC106") ~ 6,
        str_detect(arm, "LC115") ~ 15,
        TRUE ~ NA_real_
      ),
    n_detected_among_expected = 
      case_when(
        arm == "Pl" ~ n_LBP,
        str_detect(arm, "LC106") ~ n_LC106,
        str_detect(arm, "LC115") ~ n_LBP,
        TRUE ~ NA_real_
      ),
    p_detected = n_detected_among_expected/n_expected
  )

df_summ |> 
  filter(randomized) |> 
  ggplot() +
  aes(x = visit_code, y = pid, fill = total_LBP |> log10()) +
  ggh4x::facet_nested(arm + site ~ ., scales = "free_y", space = "free_y") +
  geom_tile() +
  geom_text(aes(label = n_LBP, col = total_LBP |> log10()), size = 2) +
  scale_fill_viridis_c(option = "F") + guides(col = "none") +
  scale_color_viridis_c(option = "F", direction = -1, na.value = "white") +
  theme(
    axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
    legend.position = "bottom", legend.direction = "horizontal"
  ) +
  labs(caption = "Text = number of detected (>0) LBP strains by qPCR.")


df_summ |> 
  filter(randomized) |> 
  ggplot() +
  aes(x = visit_code, y = pid, fill = n_LBP) +
  ggh4x::facet_nested(arm + site ~ ., scales = "free_y", space = "free_y") +
  geom_tile() +
  scale_fill_viridis_c("") +
  theme(
    axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
    legend.position = "bottom", legend.direction = "horizontal"
  ) +
  labs(caption = "Number of detected (>0) LBP strains by qPCR.")


df_summ |> 
  filter(randomized) |> 
  ggplot() +
  aes(x = visit_code, y = pid, fill = p_detected) +
  ggh4x::facet_nested(arm + site ~ ., scales = "free_y", space = "free_y") +
  geom_tile() +
  # scale_fill_gradient2(low = "gray80", mid = "purple", high = "blue", midpoint = 7.5) +
  scale_fill_viridis_c("proportion of LBP strains\ndetected among expected", option = "G") +
  theme(
    axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
    legend.position = "bottom", legend.direction = "horizontal"
  ) +
  labs(caption = ".")

Comparison with self-declared exposure.

Dots indicate whether participants reported inserting the study product on each day.

df_summ <- 
  df_summ |> 
  select(-starts_with("study_product")) |> 
  left_join(
    mae@metadata$exposures |> 
      select(uid, pid, visit_code, starts_with("study_product"))
  ) 

df_summ |> 
  filter(randomized) |> 
  ggplot() +
  aes(x = visit_code, y = pid, fill = total_LBP |> log10()) +
  ggh4x::facet_nested(arm + site ~ ., scales = "free_y", space = "free_y") +
  geom_tile() +
  geom_point(aes(shape = study_product |> factor()), size = 2, fill = "black", col = "white") +
  scale_shape_manual("self-declared\nexposure", values = c(NA, 21)) +
  scale_fill_viridis_c(option = "F") + 
  theme(
    axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
    legend.position = "bottom", legend.direction = "horizontal"
  ) 

df_summ |> 
  filter(randomized) |> 
  ggplot() +
  aes(x = visit_code, y = pid, fill = p_detected) +
  ggh4x::facet_nested(arm + site ~ ., scales = "free_y", space = "free_y") +
  geom_tile() + # 
  geom_point(aes(shape = study_product |> factor()), size = 2, fill = "black", col = "white") +
  scale_shape_manual("self-declared\nexposure", values = c(NA, 21)) +
  scale_fill_viridis_c(
    "proportion of LBP strains\ndetected among expected", option = "G"
    ) +
  theme(
    axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
    legend.position = "bottom", legend.direction = "horizontal"
  )

Here we adjust for reported product losses, representing actual exposure (i.e., tablets that remained in place).

df_summ |> 
  filter(randomized) |> 
  ggplot() +
  aes(x = visit_code, y = pid, fill = total_LBP |> log10()) +
  ggh4x::facet_nested(arm + site ~ ., scales = "free_y", space = "free_y") +
  geom_tile() +
  geom_point(aes(shape = study_product_PP |> round() |> factor()), size = 2, fill = "black", col = "white") +
  scale_shape_manual("self-declared\nexposure", values = c(NA, 21)) +
  scale_fill_viridis_c(option = "F") + 
  theme(
    axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
    legend.position = "bottom", legend.direction = "horizontal"
  )

df_summ |> 
  filter(randomized) |> 
  ggplot() +
  aes(x = visit_code, y = pid, fill = p_detected) +
  ggh4x::facet_nested(arm + site ~ ., scales = "free_y", space = "free_y") +
  geom_tile() + # 
  geom_point(aes(shape = study_product_PP |> round() |> factor()), size = 2, fill = "black", col = "white") +
  scale_shape_manual("self-declared\nexposure", values = c(NA, 21)) +
  scale_fill_viridis_c(
    "proportion of LBP strains\ndetected among expected", option = "G"
    ) +
  theme(
    axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
    legend.position = "bottom", legend.direction = "horizontal"
  ) 

Comparison with applicator staining

tmp <- 
   metadata(mae)$participant_crfs_merged |> 
   as_tibble() |> 
   select(pid, used_applicators, number_applicators, applicator_stain_positive, applicator_stain_negative, applicator_stain_indeterminant, total_applicators_stained, proportion_positive_applicators, total) |> 
   dplyr::rename(
     returned_used_applicators = used_applicators,
     number_of_returned_applicators = number_applicators
   )

applicator <- 
  mae@colData |> as_tibble() |>
  select(pid, site, randomized, arm, ITT, mITT, PP) |>
  distinct() |>
  filter(randomized, mITT) |>
  left_join(tmp, by = c("pid")) 


applicator <- 
  applicator |> 
  mutate(
    returned_complete_set = (number_of_returned_applicators >= 7),
    return_category = case_when(
      (returned_used_applicators == "Yes") & returned_complete_set ~ "Returned complete set",
      (returned_used_applicators == "Yes") & !returned_complete_set ~ "Returned incomplete set",
      returned_used_applicators == "No" ~ "Did not return applicators",
      TRUE ~ "Missing"
    ) |> fct_infreq()
  ) 

applicator |> 
  dplyr::count(return_category) |> 
  mutate(`%` = round(n / sum(n), 3) * 100) |> 
  gt() |> 
  cols_label(return_category = "")

	n	%
Returned complete set	77	85.6
Returned incomplete set	7	7.8
Did not return applicators	6	6.7

applicator <- 
  applicator |> 
  mutate(
    staining_done_category = 
      case_when(
        (number_of_returned_applicators == 0) | is.na(number_of_returned_applicators) ~ "No applicators returned",
        is.na(total_applicators_stained) ~ "No staining done",
        total_applicators_stained == number_of_returned_applicators ~ "All returned applicators stained",
        total_applicators_stained < number_of_returned_applicators ~ "Partial staining",
        total_applicators_stained > number_of_returned_applicators ~ "???",
        TRUE ~ "???"
      )
  )

df_summ <- 
  df_summ |> 
  select(-contains("applicator")) |> 
  left_join(
    applicator |> 
      select(
        pid, returned_used_applicators, number_of_returned_applicators, 
             starts_with("applicator_stain")
        )
  ) |> 
  mutate(
    number_of_returned_applicators = number_of_returned_applicators |> replace_na(0),
    n_pos_applicator = applicator_stain_positive |> replace_na(0),
    n_undetermined_applicator = applicator_stain_indeterminant |> replace_na(0),
    p_pos_applicator = ifelse(number_of_returned_applicators == 0, 0.4, n_pos_applicator/number_of_returned_applicators),
    p_pos_or_undet_applicator = 
      ifelse(
        number_of_returned_applicators == 0, 0.4, 
        (n_pos_applicator + n_undetermined_applicator)/number_of_returned_applicators
      )
  )

df_summ |>
  filter(randomized) |> 
  ggplot() +
  aes(x = visit_code, y = pid, fill = total_LBP |> log10()) +
  ggh4x::facet_nested(arm + site ~ ., scales = "free_y", space = "free_y") +
  geom_tile() +
  scale_fill_viridis_c(option = "F") + 
  ggnewscale::new_scale_fill() +
  geom_point(
    data = df_summ |> filter(randomized, study_product == 1),
    aes(shape = returned_used_applicators, fill = p_pos_applicator), 
    size = 2, col = "white"
    ) +
  scale_fill_gradient2(
    "Prop. pos.\napplicators",
    high = "black", mid = "green", low = "yellow", 
    midpoint = 0.5, breaks = c(0, 1/2, 1)
    )  +
  scale_shape_manual(
    "Returned\napplicators", values = c(23, 21), 
    guide = guide_legend(override.aes = list(color = "black", fill = "black"))
    ) +
  theme(
    axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
    legend.position = "bottom", legend.direction = "horizontal"
  ) 


df_summ |> 
  filter(randomized) |> 
  ggplot() +
  aes(x = visit_code, y = pid, fill = p_detected) +
  ggh4x::facet_nested(arm + site ~ ., scales = "free_y", space = "free_y") +
  geom_tile() + 
  scale_fill_viridis_c(
    "proportion of LBP strains\ndetected among expected", option = "G"
    ) +
  ggnewscale::new_scale_fill() +
  geom_point(
    data = df_summ |> filter(randomized, study_product == 1),
    aes(shape = returned_used_applicators, fill = p_pos_applicator), 
    size = 2, col = "white"
  ) +
  scale_fill_gradient2(
    "Prop. pos.\napplicators",
    high = "black", mid = "green", low = "yellow", 
    midpoint = 0.5, breaks = c(0, 1/2, 1)
    ) +
  scale_shape_manual(
    "Returned\napplicators", values = c(23, 21), 
    guide = guide_legend(override.aes = list(color = "black", fill = "black"))
    ) +
  theme(
    axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
    legend.position = "bottom", legend.direction = "horizontal"
  ) 

df_summ |> 
  filter(randomized) |> 
  ggplot() +
  aes(x = visit_code, y = pid, fill = total_LBP |> log10()) +
  ggh4x::facet_nested(arm + site ~ ., scales = "free_y", space = "free_y") +
  geom_tile() +
  scale_fill_viridis_c(option = "F") + 
  ggnewscale::new_scale_fill() +
  geom_point(
    data = df_summ |> filter(randomized, study_product == 1),
    aes(shape = returned_used_applicators, fill = p_pos_or_undet_applicator), 
    size = 2, col = "white"
    ) +
  scale_fill_gradient2(
    "Prop. pos. or ?\napplicators",
    high = "black", mid = "green", low = "yellow", 
    midpoint = 0.5, breaks = c(0, 1/2, 1)
    )  +
  scale_shape_manual(
    "Returned\napplicators", values = c(23, 21), 
    guide = guide_legend(override.aes = list(color = "black", fill = "black"))
    ) +
  theme(
    axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
    legend.position = "bottom", legend.direction = "horizontal"
  ) 

df_summ |> 
  filter(randomized) |> 
  ggplot() +
  aes(x = visit_code, y = pid, fill = p_detected) +
  ggh4x::facet_nested(arm + site ~ ., scales = "free_y", space = "free_y") +
  geom_tile() + 
  scale_fill_viridis_c(
    "proportion of LBP strains\ndetected among expected", option = "G"
    ) +
  ggnewscale::new_scale_fill() +
  geom_point(
    data = df_summ |> filter(randomized, study_product == 1),
    aes(shape = returned_used_applicators, fill = p_pos_or_undet_applicator), 
    size = 2, col = "white"
  ) +
  scale_fill_gradient2(
    "Prop. pos. or ?\napplicators",
    high = "black", mid = "green", low = "yellow", 
    midpoint = 0.5, breaks = c(0, 1/2, 1)
    ) +
  scale_shape_manual(
    "Returned\napplicators", values = c(23, 21), 
    guide = guide_legend(override.aes = list(color = "black", fill = "black"))
    ) +
  theme(
    axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
    legend.position = "bottom", legend.direction = "horizontal"
  ) 

Exposed by qPCR

For each study day, we define “effective exposure by qPCR” as having at least 50% of expected LBP strains detected by qPCR and a total LBP load exceeding 10⁷ copies per swab.

df_summ <- 
  df_summ |> 
  mutate(
    rel_ab = pmin(total_LBP/total_16S/2, 1),
    exposed = (p_detected >= 0.5) & (total_LBP > 10^7)
  )

df_summ |> 
  filter(randomized) |> 
  ggplot() +
  aes(x = visit_code, y = pid, fill = total_LBP |> log10()) +
  ggh4x::facet_nested(arm + site ~ ., scales = "free_y", space = "free_y") +
  geom_tile() +
  geom_text(aes(label = n_LBP, col = total_LBP |> log10()), size = 2) +
  scale_fill_viridis_c(option = "F") + guides(col = "none") +
  scale_color_viridis_c(option = "F", direction = -1, na.value = "white") +
  theme(
    axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
    legend.position = "bottom", legend.direction = "horizontal"
  ) +
  labs(caption = "Text = number of detected (>0) LBP strains by qPCR.")


df_summ |> 
  filter(randomized) |> 
  ggplot() +
  aes(x = visit_code, y = pid, fill = rel_ab) +
  ggh4x::facet_nested(arm + site ~ ., scales = "free_y", space = "free_y") +
  geom_tile() +
  scale_fill_viridis_c(option = "F") + guides(col = "none") +
  theme(
    axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
    legend.position = "bottom", legend.direction = "horizontal"
  ) +
  labs(caption = "Text = number of detected (>0) LBP strains by qPCR.")



df_summ |> 
  filter(randomized) |> 
  ggplot() +
  aes(x = visit_code, y = pid, fill = n_LBP) +
  ggh4x::facet_nested(arm + site ~ ., scales = "free_y", space = "free_y") +
  geom_tile() +
  scale_fill_viridis_c("") +
  theme(
    axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
    legend.position = "bottom", legend.direction = "horizontal"
  ) +
  labs(caption = "Number of detected (>0) LBP strains by qPCR.")


df_summ |> 
  filter(randomized) |> 
  ggplot() +
  aes(x = visit_code, y = pid, fill = p_detected) +
  ggh4x::facet_nested(arm + site ~ ., scales = "free_y", space = "free_y") +
  geom_tile() +
  # scale_fill_gradient2(low = "gray80", mid = "purple", high = "blue", midpoint = 7.5) +
  scale_fill_viridis_c("proportion of LBP strains\ndetected among expected", option = "G") +
  theme(
    axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
    legend.position = "bottom", legend.direction = "horizontal"
  ) +
  labs(caption = ".")

df_summ |> 
  filter(randomized) |> 
  ggplot() +
  aes(x = visit_code, y = pid, fill = exposed) +
  ggh4x::facet_nested(arm + site ~ ., scales = "free_y", space = "free_y") +
  geom_tile() +
  scale_fill_manual(
    "Exposed by qPCR",
    values = c("TRUE" = "blue", "FALSE" = "gray80"),
    labels = c("TRUE" = "Yes", "FALSE" = "No")
    ) +
  theme(
    axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
    legend.position = "bottom", legend.direction = "horizontal"
  ) +
  labs(caption = ".")

Finally, we summarize, for each participant, the number and proportion of days exposed by qPCR (daily home swabs) among expected dosing days

df_summ_pid <- 
  df_summ |> 
  filter(randomized) |> 
  mutate(
    dosing_days = 
      case_when(
        arm %in% c("Pl","LC106-7", "LC115") ~ visit_code %in% c(1102:1201), # "LC106-3", 
        arm == "LC106-3" ~ visit_code %in% c(1102:1104), 
        arm == "LC106-o" ~ visit_code %in% c(1004:1103) 
      ),
    exposed_tmp = 
      case_when(
        arm %in% c("Pl","LC106-7", "LC106-3", "LC115") ~ exposed & dosing_days,
        arm == "LC106-o" ~ exposed & dosing_days
      )
  ) |> 
  group_by(arm, site, pid, ITT, mITT, PP) |> 
  summarise(
    pid_exposed = any(exposed_tmp, na.rm = TRUE),
    n_days_exposed = sum(exposed_tmp, na.rm = TRUE),
    n_days_dosed = sum(dosing_days, na.rm = TRUE),
    .groups = "drop"
  ) |> 
  mutate(
    prop_days_exposed = ifelse(n_days_exposed == 0, 0, n_days_exposed/n_days_dosed)
  )

df_summ_pid <- 
  df_summ_pid |> 
  mutate(
    pid_exposed_by_qPCR = pid_exposed,
    n_days_exposed_by_qPCR = n_days_exposed,
    n_days_exp_dosed = n_days_dosed
  )

We add these variables to the calculated variables table in the MAE metadata

selected_vars <- c("pid_exposed_by_qPCR", "n_days_exposed_by_qPCR", "n_days_exp_dosed")

mae@metadata$calc_vars_pid <- 
  mae@metadata$calc_vars_pid |> 
  select(-any_of(selected_vars)) |> 
  left_join(
    df_summ_pid |> 
      select(pid, all_of(selected_vars)),
    by = join_by(pid)
  )

Saving the MultiAssayExperiment objects

The unblinded and augmented MAE object is saved to disk as an .rds file that can be loaded for downstream analyses.

saveRDS(
  mae, 
  str_c(
    get_05_output_dir(),  
    "mae_full_", today() |> str_remove_all("-"), ".rds"
  )
)